function [C_f f]=Push_Relabel(C);
%[C_f f]=Push_Relabel(C);
%C:capacity
%C_f:residual network
%f:maxflow

C_size=size(C);
if C_size(1)==0 || C_size(1)~=C_size(2)
    C_f=[];
    f=[];
    return;
end
V=C_size(1);
V_data=zeros(V,2);
%index define
H=1;%height
E=2;%esceed
%init V_data
C_f=C;
V_data(1,H)=V;
V_data(1,E)=-sum(C_f(1,:),2);
V_data(2:end,E)=C_f(1,2:end)';
C_f(1,2:end)=C_f(1,2:end)-V_data(2:end,E)';
C_f(2:end,1)=C_f(2:end,1)+V_data(2:end,E);
while 1
    %find a vertex which is overflow
    V_idx=find(V_data(1:end-1,E)>0,1);
    if isempty(V_idx)
        break;
    end
    %push or relabel
    V_accessible=find(C_f(V_idx,:)>0);
    [min_H min_idx]=min(V_data(V_accessible,H));
    if min_H<V_data(V_idx,H)%push
        push_idx=V_accessible(min_idx);
        push_value=min(V_data(V_idx,E),C_f(V_idx,push_idx));
        V_data(push_idx,E)=V_data(push_idx,E)+push_value;
        V_data(V_idx,E)=V_data(V_idx,E)-push_value;
        C_f(V_idx,push_idx)=C_f(V_idx,push_idx)-push_value;
        C_f(push_idx,V_idx)=C_f(push_idx,V_idx)+push_value;
    else%relabel
        V_data(V_idx,H)=min_H+1;
    end
end
f=C-C_f;
f(f<0)=0;
end